R Notebook
This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed. ## Useful Resources - Matt Ashby Crime Mapping course: https://github.com/mpjashby/crimemapping/ - Spatial Modelling for Data Scientists: https://gdsl-ul.github.io/san/ - R for Data Science: https://r4ds.had.co.nz/index.html - Geocomputation with R: https://geocompr.robinlovelace.net/
Exploring Predictors of COVID Crime Shifts
In this notebook, I’ll be prediting crime trends in London by MSOA, and looking at where COVID has been most impactive. I’ll then try to find correlates.
Tasks
- Read all burglary and robbery data - 2018 onwards
- Assign to MSOA
- Predict trend by MSOA
- Identify error by MSOA
- Visualise and model
# Data manipulation, transformation and visualisation
library(tidyverse)
# Nice tables
library(kableExtra)
# Simple features (a standardised way to encode vector data ie. points, lines, polygons)
library(sf)
# Spatial objects conversion
library(sp)
# Thematic maps
library(tmap)
# Colour palettes
library(RColorBrewer)
# More colour palettes
library(viridis)
library(raster) # raster data
library(rgdal) # input/output, projections
library(rgeos) # geometry ops
library(spdep) # spatial dependence
library(Metrics)
library(caret)Loading required package: lattice
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Attaching package: 㤼㸱caret㤼㸲
The following objects are masked from 㤼㸱package:Metrics㤼㸲:
precision, recall
The following object is masked from 㤼㸱package:purrr㤼㸲:
liftOne of the first things I notice is that while Python code is generally quite careful about imports, here we globally import everything…which is nice, but I’m also not quite clear which functions are coming from which libraries.
Now, let’s import all our crime data. Let’s start by one dataframe before figuring how to automate and concatenate. Notice R uses slashes that are the otehr way to Python and Windows
test_df <- read.csv("crimes/2018-01/2018-01-metropolitan-street.csv")
test_dfNALet’s look at all the unique crime types. Notice how we have a similar to the “unique” method in Python, and access a specific function using the same syntax.
unique(test_df["Crime.type"])To avoid this getting particularly computationally intensive, let’s write a function to pull out robberies and burglaries, and assign them a specific MSOA. Then we can iterate over all our months and get monthly counts for each offence type.
subset_df <- filter(test_df, Crime.type=="Burglary" | Crime.type=="Robbery")
subset_dfWeirdly, here I don’t need quotation marks for my colum name…not sure what’s driving that. But it’s easy to select a subset by column, and use logical comparators.
With that in mind, let’s now read our MSOA borders and assign these to an MSOA. I’ve set the OSGB CRS code as it doesn’t seem to automatically assign it.
lsoa_borders <- st_read("msoa_borders/MSOA_2011_London_gen_MHW.tab", crs=27700)Reading layer `MSOA_2011_London_gen_MHW' from data source `C:\Users\andre\Dropbox\Data Projects\Covid_crime_shift\msoa_borders\MSOA_2011_London_gen_MHW.tab' using driver `MapInfo File'
Simple feature collection with 983 features and 12 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Gridlsoa_bordersSimple feature collection with 983 features and 12 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Grid
First 10 features:
MSOA11CD MSOA11NM LAD11CD LAD11NM RGN11CD RGN11NM UsualRes HholdRes ComEstRes PopDen Hholds AvHholdSz
1 E02000001 City of London 001 E09000001 City of London E12000007 London 7375 7187 188 25.5 4385 1.6
2 E02000002 Barking and Dagenham 001 E09000002 Barking and Dagenham E12000007 London 6775 6724 51 31.3 2713 2.5
3 E02000003 Barking and Dagenham 002 E09000002 Barking and Dagenham E12000007 London 10045 10033 12 46.9 3834 2.6
4 E02000004 Barking and Dagenham 003 E09000002 Barking and Dagenham E12000007 London 6182 5937 245 24.8 2318 2.6
5 E02000005 Barking and Dagenham 004 E09000002 Barking and Dagenham E12000007 London 8562 8562 0 72.1 3183 2.7
6 E02000007 Barking and Dagenham 006 E09000002 Barking and Dagenham E12000007 London 8791 8672 119 50.6 3441 2.5
7 E02000008 Barking and Dagenham 007 E09000002 Barking and Dagenham E12000007 London 11569 11564 5 81.5 4591 2.5
8 E02000009 Barking and Dagenham 008 E09000002 Barking and Dagenham E12000007 London 8395 8376 19 87.4 3212 2.6
9 E02000010 Barking and Dagenham 009 E09000002 Barking and Dagenham E12000007 London 8615 8615 0 76.8 3292 2.6
10 E02000011 Barking and Dagenham 010 E09000002 Barking and Dagenham E12000007 London 6187 6086 101 38.8 2289 2.7
geometry
1 MULTIPOLYGON (((532135.1 18...
2 MULTIPOLYGON (((548881.6 19...
3 MULTIPOLYGON (((549102.4 18...
4 MULTIPOLYGON (((551550 1873...
5 MULTIPOLYGON (((549099.6 18...
6 MULTIPOLYGON (((549819.9 18...
7 MULTIPOLYGON (((548171.4 18...
8 MULTIPOLYGON (((546855 1863...
9 MULTIPOLYGON (((549618.8 18...
10 MULTIPOLYGON (((550244.1 18...Notice that when you’re reading a frame that isn’t “tidy”, it’s messy as hell - this is a geodataframe. Makes me miss the elegance of Geopandas somewhat. Still, it’s easy to import, and you get pretty plots. Given MOPAC data is in national grid, we’re also going to have to re-project this.
It’s also very easy to plot.
plot(lsoa_borders)plotting the first 9 out of 12 attributes; use max.plot = 12 to plot all
No idea why it wants that much white space though….
Now, let’s reproject and do a spatial join to assign all of my crimes to an MSOA. Given we’re working on London, let’s change change everything to that. We’ll start by changing my crime dataframe, that currently has latitude and longitudes as just numbers, to a spatial dataframe with coordinates
'subset_spatial <- st_as_sf(subset_df, coords = c("Longitude", "Latitude"),
crs = 4326, remove = FALSE)
subset_spatial'[1] "subset_spatial <- st_as_sf(subset_df, coords = c(\"Longitude\", \"Latitude\"), \n crs = 4326, remove = FALSE)\n\nsubset_spatial"Ah, we have missing values. Time to learn how to drop those.
In Pandas, we have easy functions to “drop_na” and “is_na” - I’m hoping to quickly find equivalents. My favourite Python approach to this is df.isna().sum(), counting how many “true” values you have when filtering like that. Can we duplicate that process?
sum(is.na(subset_df["Longitude"]))[1] 82sum(is.na(subset_df["Latitude"]))[1] 82We can! Glorious. We have 82 missing coordinates. Let’s drop all those rows.
clean_df <- subset_df[!rowSums(is.na(subset_df["Longitude"])), ]
clean_dfAnnoyingly, while R does have a drop_na function, it doesn’t take a “subset” argument like Python, which means this slightly painful fudge.
We should now be able to form our spatial df.
subset_spatial <- st_as_sf(clean_df, coords = c("Longitude", "Latitude"),
crs = 4326, remove = FALSE)
subset_spatialSimple feature collection with 10419 features and 12 fields
geometry type: POINT
dimension: XY
bbox: xmin: -0.492381 ymin: 51.28683 xmax: 0.273434 ymax: 51.68564
geographic CRS: WGS 84
First 10 features:
Crime.ID Month Reported.by Falls.within Longitude Latitude
1 628e0d673aa1b6a70479342a64b02884499df85b18dcd63cc9bff3cff9f704bc 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
2 f8e9db16dca534a83493198a838567aa5adc9dd56496edc2fff5bb4c62b8303e 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
3 cc34822074b130f141f16d02fdb2d500c86e22ae18324b43a3231b381af3f45c 2018-01 Metropolitan Police Service Metropolitan Police Service 0.135554 51.58499
4 10de581c3cd0a8c9b970824cd7589d13148d63a70b3115d95ef6c24dc0bd2c3b 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
5 50ad5d2dfea24afec9e17218db62b3d29786775db1060634ae7d4a6e7cafc3ff 2018-01 Metropolitan Police Service Metropolitan Police Service 0.127794 51.58419
6 95abc6eb0b755c9250d19bbe0062fcd4a509b701964d89667401c9dc96ca257d 2018-01 Metropolitan Police Service Metropolitan Police Service 0.138439 51.57850
7 035cc894d732addb5009148d8e163e6360094cfe451f621348f1c0419b9cbc77 2018-01 Metropolitan Police Service Metropolitan Police Service 0.138439 51.57850
8 495cac920dcf9e0e4927074e8ac307f17d340f01c69e434c4a3721df017cd342 2018-01 Metropolitan Police Service Metropolitan Police Service 0.139479 51.57974
9 48234e70cbc22265ee7968da92df1ca72f83b45414cce486ec2203daa4e59fa2 2018-01 Metropolitan Police Service Metropolitan Police Service 0.135119 51.57849
10 3f4ba20780987c37816ff34fd0f7760cf503b2e648777f84ee0104432cb01d66 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140452 51.58110
Location LSOA.code LSOA.name Crime.type Last.outcome.category Context geometry
1 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Offender sent to prison NA POINT (0.140035 51.58911)
2 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Investigation complete; no suspect identified NA POINT (0.140035 51.58911)
3 On or near Rose Lane E01000027 Barking and Dagenham 001A Burglary Status update unavailable NA POINT (0.135554 51.58499)
4 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Status update unavailable NA POINT (0.140035 51.58911)
5 On or near Hope Close E01000028 Barking and Dagenham 001B Burglary Status update unavailable NA POINT (0.127794 51.58419)
6 On or near Geneva Gardens E01000029 Barking and Dagenham 001C Burglary Investigation complete; no suspect identified NA POINT (0.138439 51.5785)
7 On or near Geneva Gardens E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA POINT (0.138439 51.5785)
8 On or near Yew Tree Gardens E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA POINT (0.139479 51.57974)
9 On or near Portland Close E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA POINT (0.135119 51.57849)
10 On or near Pedestrian Subway E01000030 Barking and Dagenham 001D Robbery Investigation complete; no suspect identified NA POINT (0.140452 51.5811)plot(subset_spatial)plotting the first 9 out of 12 attributes; use max.plot = 12 to plot all
Success! That looks faintly promising. Now, let’s figure out how to re-project.
latlong = "+init=epsg:4326"
ukgrid = "+init=epsg:27700"
subset_osgb <- st_transform(subset_spatial, ukgrid)GDAL Message 1: +init=epsg:XXXX syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order.subset_osgbSimple feature collection with 10419 features and 12 fields
geometry type: POINT
dimension: XY
bbox: xmin: 504499 ymin: 155908 xmax: 557677 ymax: 200168
projected CRS: OSGB 1936 / British National Grid
First 10 features:
Crime.ID Month Reported.by Falls.within Longitude Latitude
1 628e0d673aa1b6a70479342a64b02884499df85b18dcd63cc9bff3cff9f704bc 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
2 f8e9db16dca534a83493198a838567aa5adc9dd56496edc2fff5bb4c62b8303e 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
3 cc34822074b130f141f16d02fdb2d500c86e22ae18324b43a3231b381af3f45c 2018-01 Metropolitan Police Service Metropolitan Police Service 0.135554 51.58499
4 10de581c3cd0a8c9b970824cd7589d13148d63a70b3115d95ef6c24dc0bd2c3b 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
5 50ad5d2dfea24afec9e17218db62b3d29786775db1060634ae7d4a6e7cafc3ff 2018-01 Metropolitan Police Service Metropolitan Police Service 0.127794 51.58419
6 95abc6eb0b755c9250d19bbe0062fcd4a509b701964d89667401c9dc96ca257d 2018-01 Metropolitan Police Service Metropolitan Police Service 0.138439 51.57850
7 035cc894d732addb5009148d8e163e6360094cfe451f621348f1c0419b9cbc77 2018-01 Metropolitan Police Service Metropolitan Police Service 0.138439 51.57850
8 495cac920dcf9e0e4927074e8ac307f17d340f01c69e434c4a3721df017cd342 2018-01 Metropolitan Police Service Metropolitan Police Service 0.139479 51.57974
9 48234e70cbc22265ee7968da92df1ca72f83b45414cce486ec2203daa4e59fa2 2018-01 Metropolitan Police Service Metropolitan Police Service 0.135119 51.57849
10 3f4ba20780987c37816ff34fd0f7760cf503b2e648777f84ee0104432cb01d66 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140452 51.58110
Location LSOA.code LSOA.name Crime.type Last.outcome.category Context geometry
1 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Offender sent to prison NA POINT (548349 189976)
2 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Investigation complete; no suspect identified NA POINT (548349 189976)
3 On or near Rose Lane E01000027 Barking and Dagenham 001A Burglary Status update unavailable NA POINT (548052 189507.9)
4 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Status update unavailable NA POINT (548349 189976)
5 On or near Hope Close E01000028 Barking and Dagenham 001B Burglary Status update unavailable NA POINT (547517 189404)
6 On or near Geneva Gardens E01000029 Barking and Dagenham 001C Burglary Investigation complete; no suspect identified NA POINT (548273 188793)
7 On or near Geneva Gardens E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA POINT (548273 188793)
8 On or near Yew Tree Gardens E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA POINT (548341 188933)
9 On or near Portland Close E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA POINT (548043 188785)
10 On or near Pedestrian Subway E01000030 Barking and Dagenham 001D Robbery Investigation complete; no suspect identified NA POINT (548404 189086)Error messages in R are definitely harder to digest for me so far…I’m hoping that will pass with time. I’m also finding the documentation slightly harder to figure out, with fewer worked examples. Still, so far, so easily translateable! Now, let’s spatial join this up.
crime_with_msoa <- st_join(subset_osgb, lsoa_borders["MSOA11CD"])
crime_with_msoaSimple feature collection with 10419 features and 13 fields
geometry type: POINT
dimension: XY
bbox: xmin: 504499 ymin: 155908 xmax: 557677 ymax: 200168
projected CRS: OSGB 1936 / British National Grid
First 10 features:
Crime.ID Month Reported.by Falls.within Longitude Latitude
1 628e0d673aa1b6a70479342a64b02884499df85b18dcd63cc9bff3cff9f704bc 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
2 f8e9db16dca534a83493198a838567aa5adc9dd56496edc2fff5bb4c62b8303e 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
3 cc34822074b130f141f16d02fdb2d500c86e22ae18324b43a3231b381af3f45c 2018-01 Metropolitan Police Service Metropolitan Police Service 0.135554 51.58499
4 10de581c3cd0a8c9b970824cd7589d13148d63a70b3115d95ef6c24dc0bd2c3b 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140035 51.58911
5 50ad5d2dfea24afec9e17218db62b3d29786775db1060634ae7d4a6e7cafc3ff 2018-01 Metropolitan Police Service Metropolitan Police Service 0.127794 51.58419
6 95abc6eb0b755c9250d19bbe0062fcd4a509b701964d89667401c9dc96ca257d 2018-01 Metropolitan Police Service Metropolitan Police Service 0.138439 51.57850
7 035cc894d732addb5009148d8e163e6360094cfe451f621348f1c0419b9cbc77 2018-01 Metropolitan Police Service Metropolitan Police Service 0.138439 51.57850
8 495cac920dcf9e0e4927074e8ac307f17d340f01c69e434c4a3721df017cd342 2018-01 Metropolitan Police Service Metropolitan Police Service 0.139479 51.57974
9 48234e70cbc22265ee7968da92df1ca72f83b45414cce486ec2203daa4e59fa2 2018-01 Metropolitan Police Service Metropolitan Police Service 0.135119 51.57849
10 3f4ba20780987c37816ff34fd0f7760cf503b2e648777f84ee0104432cb01d66 2018-01 Metropolitan Police Service Metropolitan Police Service 0.140452 51.58110
Location LSOA.code LSOA.name Crime.type Last.outcome.category Context MSOA11CD
1 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Offender sent to prison NA E02000002
2 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Investigation complete; no suspect identified NA E02000002
3 On or near Rose Lane E01000027 Barking and Dagenham 001A Burglary Status update unavailable NA E02000002
4 On or near Beansland Grove E01000027 Barking and Dagenham 001A Burglary Status update unavailable NA E02000002
5 On or near Hope Close E01000028 Barking and Dagenham 001B Burglary Status update unavailable NA E02000002
6 On or near Geneva Gardens E01000029 Barking and Dagenham 001C Burglary Investigation complete; no suspect identified NA E02000002
7 On or near Geneva Gardens E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA E02000002
8 On or near Yew Tree Gardens E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA E02000002
9 On or near Portland Close E01000029 Barking and Dagenham 001C Burglary Status update unavailable NA E02000002
10 On or near Pedestrian Subway E01000030 Barking and Dagenham 001D Robbery Investigation complete; no suspect identified NA E02000002
geometry
1 POINT (548349 189976)
2 POINT (548349 189976)
3 POINT (548052 189507.9)
4 POINT (548349 189976)
5 POINT (547517 189404)
6 POINT (548273 188793)
7 POINT (548273 188793)
8 POINT (548341 188933)
9 POINT (548043 188785)
10 POINT (548404 189086)That looks like it worked - it defaults to a left join. Now, let’s group by offence type and MSOA, so as to get a count of robbery and burglary per MSOA for this month. We should then have all the code we need to create our function.
msoa_list<- crime_with_msoa %>%
group_by(MSOA11CD, Crime.type) %>%
summarize(count_by_msoa = n())`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.msoa_listSimple feature collection with 1722 features and 3 fields
geometry type: GEOMETRY
dimension: XY
bbox: xmin: 504499 ymin: 155908 xmax: 557677 ymax: 200168
projected CRS: OSGB 1936 / British National GridIt works! We’ll need to fill every missing value with 0, drop the geometry column, and then repeat the process for every month, and then we’re in business.
We’ll have to check for any values that aren’t MSOAs that aren’t present, and if they’re not, add a 0. I’m going to do this for robbery and burglary independently.
class(msoa_list)[1] "sf" "grouped_df" "tbl_df" "tbl" "data.frame"We need to remove the geometry column. Currently, our object is a spatial dataframe (sf) and a tibble (a tidyverse specific dataframe type) which is stopping me from removing the geometry data.
msoa_pivot_tibble <- as_tibble(msoa_list)
msoa_pivot_tibbleclass(msoa_pivot_tibble)[1] "tbl_df" "tbl" "data.frame"We’ve now removed the spatial frame function, and should be able to drop the last column (geometry.)
msoa_pivot_tibble <- msoa_pivot_tibble[0:3]
msoa_pivot_tibbleWe now need to fill our missing values. Rather than iterate or similar, I’ll just add an entire df filled with 0s for both crime types, then drop any duplicates - at least, that’s how I’d do it in Python, and shall try to do here!
#creating a df with all msoa names, for robbery and burglary
msoa_zero_df_robbery <- unique(as_tibble(lsoa_borders)["MSOA11CD"])
msoa_zero_df_burglary <- unique(as_tibble(lsoa_borders)["MSOA11CD"])
#adding our crime type column
msoa_zero_df_burglary["Crime.type"] = "Burglary"
msoa_zero_df_robbery["Crime.type"] = "Robbery"
#Creating a "count" column identical to our pivot, and filling it with 0
msoa_zero_df_burglary["count_by_msoa"] = as.numeric(0)
msoa_zero_df_robbery["count_by_msoa"] = as.numeric(0)
msoa_zero_df_robberyI’m a little worried about the “dbl” class, but let’s ignore that for now. Now, we need to concatenate both, and add them to our MSOA pivot.
duplicate_concat <- rbind(msoa_zero_df_robbery, msoa_zero_df_burglary)
duplicate_concatIt seems to have worked. That said, the fact there is no clear function for concatenation (in contrast to pd.concatenate in Pandas) surprises me. Now, let’s finally concatenate everything, and remove duplicates. That should form our final monthly df, and we can then combine all our previous steps into a function.
df_with_dups <- rbind(msoa_pivot_tibble, duplicate_concat)
df_with_dupsIt’s noticeable how much harder finding documentation is for R than Pandas - while the drop_duplicates function is front and center for any searches in Python, a similar search in R reveals plenty of hacky filters, but the “distinct” function seems to be what I’m actually looking for.
#creating a filter for duplicates columns, which should ignore the first instance
dup_filters <- duplicated(df_with_dups[0:2])
monthly_df <- filter(df_with_dups, !dup_filters)
monthly_dfThat should now be all our values. As a sanity check, let’s make sure we have the right number of rows, using R’s “dim” function (equivalent to shape in Pandas) to check how many unique values we would expect.
dim(unique(as_tibble(lsoa_borders)["MSOA11CD"]))[1] * 2[1] 1966We’ve got 2 extra…a bit weird, but not end of world. Let’s leave it at that.
We now need to add our monthly date to this dataframe
#select the first unique value of months in the original dataframe
month <- unique(test_df["Month"])[1,1]
monthly_df["Month"] <- month
monthly_dfNow, let’s bring all our previous work together into a function (and fix my awkward prior msoa/lsoa typo)
#quick initial function to generate our MSOA borde spatial frame, to avoid it sitting in the initial frame and gobbling loads of memory.
generate_msoa_borders <- function(file){
msoa_borders <- st_read(file, crs=27700)
return(msoa_borders)
}
make_month_pivot <- function(file){
#define our CRS
latlong = "+init=epsg:4326"
ukgrid = "+init=epsg:27700"
#read our crime from the file
test_df <- read.csv(file)
#select only our target crime types
subset_df <- filter(test_df, Crime.type=="Burglary" | Crime.type=="Robbery")
#remove any rows with a long/lat coordinate
clean_df <- subset_df[!rowSums(is.na(subset_df["Longitude"])), ]
#generate a spatial df
subset_spatial <- st_as_sf(clean_df, coords = c("Longitude", "Latitude"),
crs = 4326, remove = FALSE)
#reproject to uk grid coords
subset_osgb <- st_transform(subset_spatial, ukgrid)
#spatially join to assign to an MSOA
crime_with_msoa <- st_join(subset_osgb, msoa_borders["MSOA11CD"])
#summarise by count of MSOA
msoa_list<- crime_with_msoa %>%
group_by(MSOA11CD, Crime.type) %>%
summarize(count_by_msoa = n())
#return to a non-geographic msoa
msoa_pivot_tibble <- as_tibble(msoa_list)
msoa_pivot_tibble <- msoa_pivot_tibble[0:3]
#creating a df with all msoa names, for robbery and burglary
msoa_zero_df_robbery <- unique(as_tibble(msoa_borders)["MSOA11CD"])
msoa_zero_df_burglary <- unique(as_tibble(msoa_borders)["MSOA11CD"])
#adding our crime type column
msoa_zero_df_burglary["Crime.type"] = "Burglary"
msoa_zero_df_robbery["Crime.type"] = "Robbery"
#Creating a "count" column identical to our pivot, and filling it with 0
msoa_zero_df_burglary["count_by_msoa"] = as.numeric(0)
msoa_zero_df_robbery["count_by_msoa"] = as.numeric(0)
duplicate_concat <- rbind(msoa_zero_df_robbery, msoa_zero_df_burglary)
df_with_dups <- rbind(msoa_pivot_tibble, duplicate_concat)
#creating a filter for duplicates columns, which should ignore the first instance
dup_filters <- duplicated(df_with_dups[0:2])
monthly_df <- filter(df_with_dups, !dup_filters)
#re-add our month column
month <- unique(test_df["Month"])[1,1]
monthly_df["Month"] <- month
return(monthly_df)
}We’ve now got our tooling for a data pipeline ready to go! We can now run this on every single month of data, and aggregate into a historical combined data-set.
Data.Police.UK comes as a bunch of nested-subdirectories…in hindsight, I probably should have looked at their API, but for now let’s power ahead and figure out how to extract a list of all the CSV files in our folder and the various sub-directories.
list.files(path = "crimes") [1] "2018-01" "2018-02" "2018-03" "2018-04" "2018-05" "2018-06" "2018-07" "2018-08" "2018-09" "2018-10" "2018-11" "2018-12" "2019-01" "2019-02" "2019-03" "2019-04"
[17] "2019-05" "2019-06" "2019-07" "2019-08" "2019-09" "2019-10" "2019-11" "2019-12" "2020-01" "2020-02" "2020-03" "2020-04" "2020-05" "2020-06" "2020-07" "2020-08"
[33] "2020-09" "2020-10" "2020-11" "2020-12"As we suspected, the nested directories cause an issue - guess we’re going to have to learn about loops in R! Let’s iterate over our list of subfolders, and re-apply the function to each.
subfolders <- list.files(path = "crimes")
file_list <- list()
for (folder in subfolders){
folder_subdir <- "crimes/"
#concatenate to get our total subfolder directory - hacky but will work here.
sub_path <- paste(folder_subdir, folder, sep="")
list.files(sub_path)
file_list <- list(file_list, paste(sub_path,"/", list.files(sub_path), sep=""))
}
The bad news is, this totally didn’t work. The good news is, it led me to the far cleaner, “recursive” version of the read files function.
list.files(path = "crimes", recursive=T) [1] "2018-01/2018-01-metropolitan-street.csv" "2018-02/2018-02-metropolitan-street.csv" "2018-03/2018-03-metropolitan-street.csv"
[4] "2018-04/2018-04-metropolitan-street.csv" "2018-05/2018-05-metropolitan-street.csv" "2018-06/2018-06-metropolitan-street.csv"
[7] "2018-07/2018-07-metropolitan-street.csv" "2018-08/2018-08-metropolitan-street.csv" "2018-09/2018-09-metropolitan-street.csv"
[10] "2018-10/2018-10-metropolitan-street.csv" "2018-11/2018-11-metropolitan-street.csv" "2018-12/2018-12-metropolitan-street.csv"
[13] "2019-01/2019-01-metropolitan-street.csv" "2019-02/2019-02-metropolitan-street.csv" "2019-03/2019-03-metropolitan-street.csv"
[16] "2019-04/2019-04-metropolitan-street.csv" "2019-05/2019-05-metropolitan-street.csv" "2019-06/2019-06-metropolitan-street.csv"
[19] "2019-07/2019-07-metropolitan-street.csv" "2019-08/2019-08-metropolitan-street.csv" "2019-09/2019-09-metropolitan-street.csv"
[22] "2019-10/2019-10-metropolitan-street.csv" "2019-11/2019-11-metropolitan-street.csv" "2019-12/2019-12-metropolitan-street.csv"
[25] "2020-01/2020-01-metropolitan-street.csv" "2020-02/2020-02-metropolitan-street.csv" "2020-03/2020-03-metropolitan-street.csv"
[28] "2020-04/2020-04-metropolitan-street.csv" "2020-05/2020-05-metropolitan-street.csv" "2020-06/2020-06-metropolitan-street.csv"
[31] "2020-07/2020-07-metropolitan-street.csv" "2020-08/2020-08-metropolitan-street.csv" "2020-09/2020-09-metropolitan-street.csv"
[34] "2020-10/2020-10-metropolitan-street.csv" "2020-11/2020-11-metropolitan-street.csv" "2020-12/2020-12-metropolitan-street.csv"First, let’s create an empty dataframe we can concatenate all teh rest to.
empty_df <- tibble(
MSOA11CD = "",
Crime.type= "",
count_by_msoa= "",
Month= ""
)
empty_dfAnnoyingly, I don’t seem to be able to cr
msoa_borders <- generate_msoa_borders("msoa_borders/MSOA_2011_London_gen_MHW.tab")Reading layer `MSOA_2011_London_gen_MHW' from data source `C:\Users\andre\Dropbox\Data Projects\Covid_crime_shift\msoa_borders\MSOA_2011_London_gen_MHW.tab' using driver `MapInfo File'
Simple feature collection with 983 features and 12 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Gridmsoa_bordersSimple feature collection with 983 features and 12 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Grid
First 10 features:
MSOA11CD MSOA11NM LAD11CD LAD11NM RGN11CD RGN11NM UsualRes HholdRes ComEstRes PopDen Hholds AvHholdSz
1 E02000001 City of London 001 E09000001 City of London E12000007 London 7375 7187 188 25.5 4385 1.6
2 E02000002 Barking and Dagenham 001 E09000002 Barking and Dagenham E12000007 London 6775 6724 51 31.3 2713 2.5
3 E02000003 Barking and Dagenham 002 E09000002 Barking and Dagenham E12000007 London 10045 10033 12 46.9 3834 2.6
4 E02000004 Barking and Dagenham 003 E09000002 Barking and Dagenham E12000007 London 6182 5937 245 24.8 2318 2.6
5 E02000005 Barking and Dagenham 004 E09000002 Barking and Dagenham E12000007 London 8562 8562 0 72.1 3183 2.7
6 E02000007 Barking and Dagenham 006 E09000002 Barking and Dagenham E12000007 London 8791 8672 119 50.6 3441 2.5
7 E02000008 Barking and Dagenham 007 E09000002 Barking and Dagenham E12000007 London 11569 11564 5 81.5 4591 2.5
8 E02000009 Barking and Dagenham 008 E09000002 Barking and Dagenham E12000007 London 8395 8376 19 87.4 3212 2.6
9 E02000010 Barking and Dagenham 009 E09000002 Barking and Dagenham E12000007 London 8615 8615 0 76.8 3292 2.6
10 E02000011 Barking and Dagenham 010 E09000002 Barking and Dagenham E12000007 London 6187 6086 101 38.8 2289 2.7
geometry
1 MULTIPOLYGON (((532135.1 18...
2 MULTIPOLYGON (((548881.6 19...
3 MULTIPOLYGON (((549102.4 18...
4 MULTIPOLYGON (((551550 1873...
5 MULTIPOLYGON (((549099.6 18...
6 MULTIPOLYGON (((549819.9 18...
7 MULTIPOLYGON (((548171.4 18...
8 MULTIPOLYGON (((546855 1863...
9 MULTIPOLYGON (((549618.8 18...
10 MULTIPOLYGON (((550244.1 18...Our MSOA border helper functions seems to work. Now, time to do the heavy lifting!
subfiles <- list.files(path = "crimes", recursive=T)
for (file in subfiles){
folder_subdir <- "crimes/"
#concatenate to get our total subfolder directory - hacky but will work here.
sub_path <- paste(folder_subdir, file, sep="")
monthly_df <- make_month_pivot(sub_path)
empty_df <- rbind(empty_df, monthly_df)
}`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.
`summarise()` has grouped output by 'MSOA11CD'. You can override using the `.groups` argument.That seems to have worked! The processing time was longer than I expected (which is probably something to do with how R stores memory) - let’s look at our previously empty dataframe.
empty_dfunique(empty_df["Month"])So we now have a combined dataframe of just under 71,000 rows, for 37 individual months beetween January 2018 and December 2020, for every robbery and burglary in London, assigned to an MSOA. I’d call that a win!
This has been somewhat more painful than I expected, so before going any further, let’s figure out how to save this file. Given it’s all strings and integers, a simple CSV should do for now.
write.csv(empty_df,"msoa_crime_matrix.csv")
Modelling - Crime Trends by MSOA
We can now move on to the fun bit - predicting the crime trend we’d expect, and then looking at how much it diverges when we hit the “pandemic disruption” period.
I’ll either be using auto-arima or Facebook’s prophet algorithm, both of which produce relatively accurate forecasts with little necessary tuning. I’m hoping 2018 through early 2020 should be sufficient to establish trends and seasonality. We can then use our error rate from March 2020 onwards as a measure of the “pandemic effect”.
As a test, let’s start by predicting one MSOA: the first in our df, “E02000001”
empty_df <- read.csv("msoa_crime_matrix.csv")
empty_df <- empty_df[2:70848,2:5]
empty_dfsingle_msoa_df <- filter(empty_df, MSOA11CD == "E02000001" & Crime.type=="Burglary")
single_msoa_dfAs we’d expect, 36 rows for 36 months. Let’s covert those rows to a date, and start making predictions.
single_msoa_df$DateString <- paste(single_msoa_df$Month, "-01")
single_msoa_dfNAConverting these to dates was harder than I’d anticipated, but the Tidyverse ecosystem does have some nifty tools!
library(lubridate)
Attaching package: 㤼㸱lubridate㤼㸲
The following objects are masked from 㤼㸱package:rgeos㤼㸲:
intersect, setdiff, union
The following objects are masked from 㤼㸱package:raster㤼㸲:
intersect, union
The following objects are masked from 㤼㸱package:base㤼㸲:
date, intersect, setdiff, unionsingle_msoa_df$DateClean <- ymd(single_msoa_df$DateString)
single_msoa_dfLet’s now build our “pre-pandemic” training set, up to February 2020, and use prophet to make some predictions.
training_set <- filter(single_msoa_df, DateClean < "2020-03-01")
training_setNAtraining_df <- tibble(
ds=training_set$DateClean,
y=training_set$count_by_msoa
)
training_dfWe can now instantiate our prophet model, and start making predictions.
library(prophet)Loading required package: Rcpp
Loading required package: rlang
Attaching package: 㤼㸱rlang㤼㸲
The following object is masked from 㤼㸱package:Metrics㤼㸲:
ll
The following objects are masked from 㤼㸱package:purrr㤼㸲:
%@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl, flatten_raw, invoke, list_along, modify, prepend, splicem <- prophet(training_df)Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.
Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.
n.changepoints greater than number of observations. Using 19We’ll predict for a period of around 3 months and then compare to what actually happened.
future <- make_future_dataframe(m, periods = 6, freq = 'month')
tail(future)Let’s get forecasting
forecast <- predict(m, future)
tail(forecast[c('ds', 'yhat', 'yhat_lower', 'yhat_upper')])NA# R
plot(m, forecast)
♠# RError: unexpected input in "\"These predictions obviously look a little silly, but the yearly trend (which is what we really wanted to get out of this) doesn’t look mad to me. I’m hoping that using all our MSOA in aggregate, we’ll get meaningful data. Firstly, we need to get our error rate.
head(forecast)
forecast$Month <- month(forecast$ds)
forecast$Year <- year(forecast$ds)
forecast
To do next.
Group prediction by month We’ll look at April and May, which were “peak” London COVID effect
forecastthis_year <- filter(forecast, Year > 2019)
peak_pandemic <- filter(this_year, Month== 4 | Month== 5 )
peak_pandemicpredictionPivot <- peak_pandemic %>%
group_by(Month) %>%
summarize(predicted_burglary = mean(yhat))
predictionPivot
Now, let’s compare that to our ACTUAL data, and get an error rate.
single_msoa_df$MonthNum <- month(single_msoa_df$DateClean)
single_msoa_df$YearNum <- year(single_msoa_df$DateClean)
this_year_actual <- filter(single_msoa_df, YearNum > 2019)
peak_pandemic_actual <- filter(this_year_actual, MonthNum== 4 | MonthNum== 5 )
peak_pandemic_actualAs such, our model anticipated around 18 burglaries during this period, while in reality, there was 1. This is the measure of our “covid error” for the MSOA. Let’s get that measure into a percentage error as well as absolute error, and then we can repeat the process for all of London
actual_burglary <- sum(peak_pandemic_actual$count_by_msoa)
pred_burglary <- sum(predictionPivot$predicted_burglary)
error <- actual_burglary - pred_burglary
percentage_error <- error / pred_burglary
print("Burglary Count")
print(actual_burglary)
print("Predicted")
print(pred_burglary)
print("Actual Error")
print(error)
print("Percentage Error")
print(percentage_error)Now, let’s automate. Similar to our previous process, we’ll create a dataframe for every MSOA, and both types, then run a for loop repeating our process across London.
msoa_error_tibble <- tibble(
MSOA11CD = "",
burglaryActual= "",
burglaryPredicted= "",
burglaryError= "",
burglaryPercentError="",
robberyActual= "",
robberyPredicted= "",
robberyError= "",
robberyPercentError=""
)
msoa_error_tibblecalculate_error <- function(msoaName){
#select only burglary and our msoa
single_msoa_df <- filter(empty_df, MSOA11CD == msoaName & Crime.type=="Burglary")
#clean date date
single_msoa_df$DateString <- paste(single_msoa_df$Month, "-01")
single_msoa_df$DateClean <- ymd(single_msoa_df$DateString)
#generate training set up until March
training_set <- filter(single_msoa_df, DateClean < "2020-03-01")
#prepare for Prophet
training_df <- tibble(
ds=training_set$DateClean,
y=training_set$count_by_msoa)
#start and predict prophet for 6 months
m <- prophet(training_df)
future <- make_future_dataframe(m, periods = 6, freq = 'month')
forecast <- predict(m, future)
forecast$Month <- month(forecast$ds)
forecast$Year <- year(forecast$ds)
#aggregate forecasts and actual crime
this_year <- filter(forecast, Year > 2019)
peak_pandemic <- filter(this_year, Month== 4 | Month== 5 )
predictionPivot <- peak_pandemic %>%
group_by(Month) %>%
summarize(predicted_burglary = mean(yhat))
single_msoa_df$MonthNum <- month(single_msoa_df$DateClean)
single_msoa_df$YearNum <- year(single_msoa_df$DateClean)
#generate error rates
this_year_actual <- filter(single_msoa_df, YearNum > 2019)
peak_pandemic_actual <- filter(this_year_actual, MonthNum== 4 | MonthNum== 5 )
actual_burglary <- sum(peak_pandemic_actual$count_by_msoa)
pred_burglary <- sum(predictionPivot$predicted_burglary)
error_burg <- actual_burglary - pred_burglary
percentage_error_burg <- error_burg / pred_burglary
#now repeat for robbery
single_msoa_df <- filter(empty_df, MSOA11CD == msoaName & Crime.type=="Robbery")
single_msoa_df$DateString <- paste(single_msoa_df$Month, "-01")
single_msoa_df$DateClean <- ymd(single_msoa_df$DateString)
training_set <- filter(single_msoa_df, DateClean < "2020-03-01")
training_df <- tibble(
ds=training_set$DateClean,
y=training_set$count_by_msoa)
m <- prophet(training_df)
future <- make_future_dataframe(m, periods = 6, freq = 'month')
forecast <- predict(m, future)
forecast$Month <- month(forecast$ds)
forecast$Year <- year(forecast$ds)
this_year <- filter(forecast, Year > 2019)
peak_pandemic <- filter(this_year, Month== 4 | Month== 5 )
predictionPivot <- peak_pandemic %>%
group_by(Month) %>%
summarize(predicted_burglary = mean(yhat))
single_msoa_df$MonthNum <- month(single_msoa_df$DateClean)
single_msoa_df$YearNum <- year(single_msoa_df$DateClean)
this_year_actual <- filter(single_msoa_df, YearNum > 2019)
peak_pandemic_actual <- filter(this_year_actual, MonthNum== 4 | MonthNum== 5 )
actual_robbery <- sum(peak_pandemic_actual$count_by_msoa)
pred_robbery <- sum(predictionPivot$predicted_burglary)
error_rob <- actual_robbery - pred_robbery
percentage_error_rob <- error_rob / pred_robbery
#create our output dataframe and return it
msoa_error_tibble <- tibble(
MSOA11CD = msoaName,
burglaryActual= actual_burglary,
burglaryPredicted= pred_burglary,
burglaryError= error_burg,
burglaryPercentError = percentage_error_burg,
robberyActual= actual_robbery,
robberyPredicted= pred_robbery,
robberyError= error_rob,
robberyPercentError=percentage_error_rob
)
return(msoa_error_tibble)
}Messy and hacky, but theoretically functional! Now for the long bit - let’s loop over all our MSOAs, and get our aggregate error dataframe.
for (msoa in unique(empty_df$MSOA11CD)){
iterated_msoa_df <- calculate_error(msoa)
msoa_error_tibble <- rbind(msoa_error_tibble, iterated_msoa_df)
}
msoa_error_tibblewrite_csv(msoa_error_tibble, "msoa_error_table.csv")msoa_error_tibblemsoa_error_tibble[,2:9] <- lapply(msoa_error_tibble[,2:9], as.numeric)msoa_error_tibble <- msoa_error_tibble[2:980, ]
msoa_error_tibbleLet’s also calculate the “Relative Percentage Difference” (RPD) of our estimates.
msoa_error_tibble$RPDBurglary <- 2*((msoa_error_tibble$burglaryPredicted - msoa_error_tibble$burglaryActual)/(abs(msoa_error_tibble$burglaryPredicted) + abs(msoa_error_tibble$burglaryActual)))
msoa_error_tibble$RPDRobbery <- 2*((msoa_error_tibble$robberyPredicted - msoa_error_tibble$robberyActual)/(abs(msoa_error_tibble$robberyPredicted) + abs(msoa_error_tibble$robberyActual)))Now, let’s link all these back to our original geographic dataframe.
lsoa_borders <- st_read("msoa_borders/MSOA_2011_London_gen_MHW.tab", crs=27700)Reading layer `MSOA_2011_London_gen_MHW' from data source `C:\Users\andre\Dropbox\Data Projects\Covid_crime_shift\msoa_borders\MSOA_2011_London_gen_MHW.tab' using driver `MapInfo File'
Simple feature collection with 983 features and 12 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Gridlsoa_bordersSimple feature collection with 983 features and 12 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Grid
First 10 features:
MSOA11CD MSOA11NM LAD11CD LAD11NM RGN11CD RGN11NM UsualRes HholdRes ComEstRes PopDen Hholds AvHholdSz
1 E02000001 City of London 001 E09000001 City of London E12000007 London 7375 7187 188 25.5 4385 1.6
2 E02000002 Barking and Dagenham 001 E09000002 Barking and Dagenham E12000007 London 6775 6724 51 31.3 2713 2.5
3 E02000003 Barking and Dagenham 002 E09000002 Barking and Dagenham E12000007 London 10045 10033 12 46.9 3834 2.6
4 E02000004 Barking and Dagenham 003 E09000002 Barking and Dagenham E12000007 London 6182 5937 245 24.8 2318 2.6
5 E02000005 Barking and Dagenham 004 E09000002 Barking and Dagenham E12000007 London 8562 8562 0 72.1 3183 2.7
6 E02000007 Barking and Dagenham 006 E09000002 Barking and Dagenham E12000007 London 8791 8672 119 50.6 3441 2.5
7 E02000008 Barking and Dagenham 007 E09000002 Barking and Dagenham E12000007 London 11569 11564 5 81.5 4591 2.5
8 E02000009 Barking and Dagenham 008 E09000002 Barking and Dagenham E12000007 London 8395 8376 19 87.4 3212 2.6
9 E02000010 Barking and Dagenham 009 E09000002 Barking and Dagenham E12000007 London 8615 8615 0 76.8 3292 2.6
10 E02000011 Barking and Dagenham 010 E09000002 Barking and Dagenham E12000007 London 6187 6086 101 38.8 2289 2.7
geometry
1 MULTIPOLYGON (((532135.1 18...
2 MULTIPOLYGON (((548881.6 19...
3 MULTIPOLYGON (((549102.4 18...
4 MULTIPOLYGON (((551550 1873...
5 MULTIPOLYGON (((549099.6 18...
6 MULTIPOLYGON (((549819.9 18...
7 MULTIPOLYGON (((548171.4 18...
8 MULTIPOLYGON (((546855 1863...
9 MULTIPOLYGON (((549618.8 18...
10 MULTIPOLYGON (((550244.1 18...geographic_error_map <- left_join(lsoa_borders, msoa_error_tibble, by = "MSOA11CD")
geographic_error_mapSimple feature collection with 983 features and 22 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 503574.2 ymin: 155850.8 xmax: 561956.7 ymax: 200933.6
projected CRS: OSGB 1936 / British National Grid
First 10 features:
MSOA11CD MSOA11NM LAD11CD LAD11NM RGN11CD RGN11NM UsualRes HholdRes ComEstRes PopDen Hholds AvHholdSz burglaryActual
1 E02000001 City of London 001 E09000001 City of London E12000007 London 7375 7187 188 25.5 4385 1.6 1
2 E02000002 Barking and Dagenham 001 E09000002 Barking and Dagenham E12000007 London 6775 6724 51 31.3 2713 2.5 8
3 E02000003 Barking and Dagenham 002 E09000002 Barking and Dagenham E12000007 London 10045 10033 12 46.9 3834 2.6 11
4 E02000004 Barking and Dagenham 003 E09000002 Barking and Dagenham E12000007 London 6182 5937 245 24.8 2318 2.6 2
5 E02000005 Barking and Dagenham 004 E09000002 Barking and Dagenham E12000007 London 8562 8562 0 72.1 3183 2.7 4
6 E02000007 Barking and Dagenham 006 E09000002 Barking and Dagenham E12000007 London 8791 8672 119 50.6 3441 2.5 5
7 E02000008 Barking and Dagenham 007 E09000002 Barking and Dagenham E12000007 London 11569 11564 5 81.5 4591 2.5 10
8 E02000009 Barking and Dagenham 008 E09000002 Barking and Dagenham E12000007 London 8395 8376 19 87.4 3212 2.6 16
9 E02000010 Barking and Dagenham 009 E09000002 Barking and Dagenham E12000007 London 8615 8615 0 76.8 3292 2.6 3
10 E02000011 Barking and Dagenham 010 E09000002 Barking and Dagenham E12000007 London 6187 6086 101 38.8 2289 2.7 2
burglaryPredicted burglaryError burglaryPercentError robberyActual robberyPredicted robberyError robberyPercentError RPDBurglary RPDRobbery
1 7.62501855 -6.6250185 -0.8688528 1 -1.976662 2.976662 -1.5059034 1.5362329 -2.0000000
2 -9.23326714 17.2332671 -1.8664322 0 1.129576 -1.129576 -1.0000000 -2.0000000 2.0000000
3 12.34800064 -1.3480006 -0.1091675 10 6.922428 3.077572 0.4445798 0.1154703 -0.3637269
4 -4.71960263 6.7196026 -1.4237645 0 -1.271292 1.271292 -1.0000000 -2.0000000 -2.0000000
5 4.58490183 -0.5849018 -0.1275713 1 9.214027 -8.214027 -0.8914698 0.1362629 1.6083817
6 13.29091514 -8.2909151 -0.6238032 5 1.267089 3.732911 2.9460534 0.9065610 -1.1912744
7 11.22961111 -1.2296111 -0.1094972 1 -5.357023 6.357023 -1.1866708 0.1158393 -2.0000000
8 -1.51347706 17.5134771 -11.5716832 3 -3.718159 6.718159 -1.8068510 -2.0000000 -2.0000000
9 -0.07818478 3.0781848 -39.3706412 0 -5.041581 5.041581 -1.0000000 -2.0000000 -2.0000000
10 -1.54900357 3.5490036 -2.2911526 1 4.929491 -3.929491 -0.7971393 -2.0000000 1.3254058
geometry
1 MULTIPOLYGON (((532135.1 18...
2 MULTIPOLYGON (((548881.6 19...
3 MULTIPOLYGON (((549102.4 18...
4 MULTIPOLYGON (((551550 1873...
5 MULTIPOLYGON (((549099.6 18...
6 MULTIPOLYGON (((549819.9 18...
7 MULTIPOLYGON (((548171.4 18...
8 MULTIPOLYGON (((546855 1863...
9 MULTIPOLYGON (((549618.8 18...
10 MULTIPOLYGON (((550244.1 18...Let’s map both of these metrics, and see what it looks like.
# map
burg_map <- tm_shape(geographic_error_map) +
tm_fill(col = "RPDRobbery", title = "Robbery Relative Error")
rob_map <-tm_shape(geographic_error_map) +
tm_fill(col = "RPDBurglary", title = "Burglary Relative Error")
tmap_arrange(burg_map, rob_map)
NA
NAAs a final part of this project, I’m going to explore some geographic modelling. Let’s start with linking our current data with the London MOPAC MSOA Atlas, which should provide a whole bunch of useful demographic and economic data. I’ve slightly modified it in Excel to get rid of the weird header structure.
library(readxl)
msoa_atlas <- read_excel("msoa_atlas/msoa-data.xls")New names:
* `House Prices Sales 2011` -> `House Prices Sales 2011...129`
* `House Prices Sales 2011` -> `House Prices Sales 2011...130`msoa_atlasLet’s do one last spatial join to bring all these things together
geographic_msoa_matrix <- left_join(geographic_error_map, msoa_atlas, by = "MSOA11CD")
Let’s provide a tible version as well for clear analysis
msoa_matrix_tbl <- as_tibble(geographic_msoa_matrix)
msoa_matrix_tblLet’s look at how correlated our factors are
library(corrr)
Attaching package: 㤼㸱corrr㤼㸲
The following object is masked from 㤼㸱package:raster㤼㸲:
stretchcorr_df <- correlate(msoa_matrix_tbl, quiet = TRUE)Error in stats::cor(x = x, y = y, use = use, method = method) :
'x' must be numericAnnoyingly, unlike Pandas, R throws up errors here (while Python implicitly gets rid of non-numerical columns - let’s clean it up
msoa_matrix_numeric <-dplyr::select_if(msoa_matrix_tbl, is.numeric)
msoa_matrix_numeric
corr_df <- correlate(dplyr::select_if(msoa_matrix_tbl, is.numeric), quiet = TRUE)
corr_dfNA
NALet’s now look for correlates of our error rate for burglary and robbery
options(scipen = 999)
dplyr::select(corr_df[order(corr_df$RPDRobbery),] , term, RPDRobbery)NAdplyr::select(corr_df[order(corr_df$RPDBurglary),] , term, RPDBurglary)There are very few decent correlates in the robbery data - everything is a bit of a mess. That’s not the case in the burglary data however: we might be able to do some modeling here. General deprivation indicators stand out very sharply, correlated to ethnicity.
I’d normally start exploring spatial models and weights, but I think that’s a little outside the scope of this first project. Instead, let’s get straight to modeling.
log(msoa_matrix_numeric["Income Deprivation (2010) % living in income deprived households reliant on means tested benefit"])msoa_burglary_copy <- msoa_matrix_numeric
names(msoa_burglary_copy)[names(msoa_burglary_copy) == "Income Deprivation (2010) % living in income deprived households reliant on means tested benefit"] <- "hhPercentBenefit"
names(msoa_burglary_copy)[names(msoa_burglary_copy) == "Lone Parents (2011 Census) Lone parents not in employment"] <- "UnempLoneParents"
names(msoa_burglary_copy)[names(msoa_burglary_copy) == "Ethnic Group (2011 Census) Black/African/Caribbean/Black British (%)"] <- "blackPercent"
names(msoa_burglary_copy)[names(msoa_burglary_copy) == "Health (2011 Census) Bad health (%)"] <- "percentBadHealth"
names(msoa_burglary_copy)[names(msoa_burglary_copy) == "Tenure (2011) Owned: Owned outright"] <- "HouseOwned"
feature_df <- dplyr::select(msoa_burglary_copy, RPDBurglary, hhPercentBenefit, UnempLoneParents, blackPercent, percentBadHealth, HouseOwned)
feature_dfWe have a few NA. Let’s get rid of them
feature_df <- drop_na(feature_df, RPDBurglary)
feature_df colnames(feature_df)[1] "RPDBurglary" "hhPercentBenefit" "UnempLoneParents" "blackPercent" "percentBadHealth" "HouseOwned" colSums(is.na(feature_df)) RPDBurglary hhPercentBenefit UnempLoneParents blackPercent percentBadHealth HouseOwned
0 0 0 0 0 0 pairs(feature_df)
We are going to want to perform a log transform. Because our RPD has negative values, we’ll need to add a constant before - in thise case, 3 (so everything is positive)
feature_df$RPDBurglaryTranform <- feature_df$RPDBurglary + 3
feature_dfNow, let’s produce our log transform columns
for (col in colnames(feature_df)){
new_name <- paste("log_", col, sep = "")
feature_df[new_name] <- log(feature_df[col])
}NaNs producedfeat_transform_df <- feature_df[,9:14]
feat_transform_dfpairs(feat_transform_df)
correlate(feat_transform_df)
Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'mod_burglary <- lm(log_RPDBurglaryTranform ~ log_blackPercent + log_hhPercentBenefit , data = feat_transform_df)
summary(mod_burglary)
Call:
lm(formula = log_RPDBurglaryTranform ~ log_blackPercent + log_hhPercentBenefit,
data = feat_transform_df)
Residuals:
Min 1Q Median 3Q Max
-1.2291 -0.0747 0.1930 0.3203 0.5826
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.26338 0.08107 15.584 <2e-16 ***
log_blackPercent -0.04064 0.03000 -1.355 0.176
log_hhPercentBenefit -0.03171 0.04473 -0.709 0.479
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4874 on 976 degrees of freedom
Multiple R-squared: 0.0121, Adjusted R-squared: 0.01008
F-statistic: 5.978 on 2 and 976 DF, p-value: 0.002627mod_burglary <- lm(log_RPDBurglaryTranform ~ log_blackPercent , data = feat_transform_df)
summary(mod_burglary)
Call:
lm(formula = log_RPDBurglaryTranform ~ log_blackPercent, data = feat_transform_df)
Residuals:
Min 1Q Median 3Q Max
-1.22673 -0.07921 0.19246 0.32606 0.56948
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.21385 0.04113 29.515 < 2e-16 ***
log_blackPercent -0.05809 0.01716 -3.385 0.00074 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4872 on 977 degrees of freedom
Multiple R-squared: 0.01159, Adjusted R-squared: 0.01058
F-statistic: 11.46 on 1 and 977 DF, p-value: 0.0007397So, we have significance. It looks like the main driver are general deprivation indicators, but due to correlation, we’re struggling to get any more detailed. Let’s use Random Forests to try and dig into this as bit further.
msoa_matrix_numericOnce again, let’s get rid of NA. This time, let’s do our RPD, and then any columns
rf_msoa_matrix <- drop_na(msoa_matrix_numeric, RPDBurglary)
rf_msoa_matrixclean_rf_matrix <- rf_msoa_matrix[ , colSums(is.na(rf_msoa_matrix)) == 0]
clean_rf_matrix